Freely decaying turbulence in two-dimensional electrostatic gyrokinetics 
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In magnetized plasmas, a turbulent cascade occurs in phase space at scales smaller than the 
thermal Larmor radius ("sub-Larmor scales") [Phys. Rev. Lett. 103, 015003 (2009)]. When the 
turbulence is restricted to two spatial dimensions perpendicular to the background magnetic field, 
two independent cascades may take place simultaneously because of the presence of two collisionless 
invariants. In the present work, freely decaying turbulence of two-dimensional electrostatic gyroki- 
netics is investigated by means of phenomenological theory and direct numerical simulations. A dual 
cascade (forward and inverse cascades) is observed in velocity space as well as in position space, 
which we diagnose by means of nonlinear transfer functions for the collisionless invariants. We find 
that the turbulence tends to a time-asymptotic state, dominated by a single scale that grows in time. 
A theory of this asymptotic state is derived in the form of decay laws. Each case that we study 
falls into one of three regimes (weakly collisional, marginal, and strongly collisional) , determined 
by a dimensionless number D*, a quantity analogous to the Reynolds number. The marginal state 
is marked by a critical number — Do that is preserved in time. Turbulence initialized above 
this value become increasingly inertial in time, evolving toward larger and larger D*; turbulence 
initialized below Do become more and more collisional, decaying to progressively smaller D*. 

PACS numbers: 52.30.Gz, 52.35.Ra, 52.65.Tt 



I. INTRODUCTION 

Plasma turbulence plays important roles in fusion de- 
vices and various space and astrophysical situations, 
where it is an essential phenomenon underlying trans- 
port of mean quantities and particle heating [ll-[ll|. For 
these collisionless or weakly collisional plasmas, such tur- 
bulence requires a kinetic description in phase space, es- 
pecially at small scales where dissipation takes place. 

Turbulence theory in kinetic phase space is more than a 
simple extension of Navier-Stokes turbulence into higher 
dimensions as velocity space is not exactly equivalent to 
position space: For instance, there is no translational 
symmetry (i.e. small velocities are not equivalent to large 
velocities), so there is always some large-scale velocity 
dependence in the problem. However, in some simpli- 
fied cases, classical fluid dynamical theories 
be naturally extended into phase space [ll|. 
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magnetized plasmas, the gyrokinetic (GK) theory plol - 
[22! provides the minimal kinetic description of the low- 
frequency turbulence. 

In electrostatic gyrokinetics, nonlinear interactions in- 
troduce a cascade of perturbed entropy (which is propor- 
tional to the perturbed free energy at the sub-Larmor 
scales) to smaller scales both in position and velocity 
space [ll|, [l^ - fTst . When the turbulence is restricted 
to two position-space dimensions, the system has two 
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collisionless invariants. One of them is the free energy 
or entropy which is also an invariant in three dimen- 
sions (3D), and another, approximately related to ki- 
netic energy in the long wave-length limit, is particular 
to two dimensions (2D). In this regard, 2D gyrokinetic 
turbulence is analogous to 2D fluid turbulence, and in- 
deed, reduces to it in a particular long wave-length limit 
[23|. These two invariants cannot share the same 
local-interaction space in a Kolmogorov-like phenomenol- 
ogy, which leads to a dual cascade (forward and inverse 
cascades) 0, [l3|, [23 - [26| . As the nonlinear term of 2D gy- 
rokinetics is identical in form to that of 3D gyrokinetics, 
the understanding of the nonlinear interaction in purely 
2D system will serve as a foundation for understanding 
general 3D magnetized plasmas [18]. 

In this paper we focus on the freely decaying turbu- 
lence problem for the electrostatic 2D GK system. We 
first introduce the GK equation briefly and describe its 
basic nature in Sec. [Ill We also review basic characteris- 
tics of the nonlinear phase mixing at sub-Larmor scales 
that are reported in [H, [13, E^- Sections [ml and |lVl 
are the main contents of the paper. In Sec. IIIIl we de- 
scribe a phenomenological theory of the dual cascade in 
freely decaying turbulence based on the theory first de- 
veloped for 2D fluid turbulence [28|. There are three 
regimes for the turbulence: These correspond to weakly 
collisional, marginal, and strongly collisional cases. In 
the strongly collisional case, dissipation acts strongly and 
the system decays to become more and more dissipative. 
In the marginal case, the collisional dissipation balances 
with nonlinear (inertial) turnover and the system is pre- 
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served in this state. For weakly collisional cases, the 
turbulence is, in a sense, able to escape the effects of dis- 
sipation, becoming less and less collisional in time, and 
tending asymptotically to a state of zero electrostatic en- 
ergy decay. Section IIVI presents the results of numerical 
simulation of the freely decaying turbulence. We demon- 
strate the inverse (forward) transfer of energy (entropy) 
in the direct numerical simulation and investigate the 
time-asymptotic decay laws of the two collisionless in- 
variants, comparing with the theory developed in Sec. 
mil We conclude with a summary of our results in Sec. 

El 



II. PHASE-SPACE TURBULENCE 

A. Gyrokinetic equations and invariants 

We first introduce the gyrokinetic (GK) model briefly 
[Hl-lill • Since we are concerned with turbulence in mag- 
netized plasmas, the dynamics of interest is much slower 
than particle gyromotion. The gyromotion is thus av- 
eraged over, eliminating gyroangle dependence from the 
system. The GK system has 3 spatial coordinates (x, 
z), and 2 velocity coordinates {y^^ v\\)^ where ± and || 
denote perpendicular and parallel directions to the back- 
ground magnetic field, respectively. We assume the back- 
ground plasma and magnetic field are uniform in space 
and time. It is necessary to distinguish between the parti- 
cle coordinate r and the gyrocenter coordinate R. These 
coordinates are connected by the Catto transform [19] 
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where is the unit vector along the background mag- 
netic field and Q is the gyrofrequency. 

We further reduce the GK equation to 2D in position 
space, or 4D in phase space [29], by ignoring variation 
along the mean field {k\\ = 0). This not only reduces 
the dimension of the system but also removes one of 
the mechanisms of creating velocity-space structure — 
linear parallel phase mixing (Landau damping), a much 
more familiar and better understood phenomenon than 
the nonlinear perpendicular phase mixing, on which we 
will concentrate in this paper. The resulting GK equa- 
tion (for the ions) is 



(2) 



where {■)r is the gyroaverage holding the guiding cen- 
ter position R constant, g = {Sf)R is the gyroaverage of 
the perturbed ion distribution function Sf , cp is the elec- 
trostatic potential, Bq is the background magnetic field 
(aligned with the z-axis), and {/, ^} = • (V/ x V^). 
The collision operator C we use in our simulations de- 
scribes pitch-angle scattering and energy diffusion with 
proper conservation properties [30|, [3l| . 



The potential (f in Eqn. (j2j) is calculated from the 
quasineutrality condition: Written in the Fourier space, 
it is 



(1 



ro)m =q^JJo mdv, (3) 



where the hat denotes the Fourier coefficients, Jo is the 
Bessel function, representing, in Fourier space, the gy- 
roaverage at fixed particle position, Fq = Io{b)e~^, Iq is 
the modified Bessel function of the first kind, b = /c^p^/2, 
p is the ion thermal Larmor radius, q is the charge, no and 
To are the density and temperature of the background 
Maxwellian Fo, and i and e are the species indices. One 
may use r = —QeTQi/qiTQe for Boltzmann-response (3D) 
electrons or r = for no-response (2D) electrons |15| . 
Hereafter we use the no-response electrons as in Refs. 
[m, |l3 because formally, the electrons cannot contribute 
to the potential if /cy = exactly [l^. This choice is not 
very important as it only introduces minor differences in 
various prefactors. 

The 2D electrostatic GK system possesses two 
quadratic positive-definite collisionless invariants [l5|. 
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There are various ways to choose two independent in- 
variants in our system. In Refs. [TgI, total perturbed 
entropy (or free energy), Wtot = W -\- E (see Eqn. (3.9) 
of Ref. p^), is used in order to make the connection to 
thermodynamics. Here we use the quantity W for the 
sake of simplicity [isj. In fact, averaged over R is 
itself conserved (that is, it is conserved for each value of 
v±) [m, [32|. However, it is sufficient for the purposes of 
this paper to consider only the integrated quantity W. 
One can adapt the arguments of Ref. [13] to show that 
the presence of conserved quantities W and E implies 
a dual cascade (i.e., both forward and inverse cascades) 

M. 



B. Nonlinear phase mixing 

Due to the neglect of the parallel streaming term, the 
creation of velocity space structure originates solely from 
the advection of the distribution function by the gyroav- 
eraged E x B drift [the nonlinear term in Eqn. (|2])]. For 
small-scale electric fields, particles with different gyro- 
radii execute different E x B motions because they "see" 
different effective potentials; this leads to nonlinear phase 
mixing and other novel phenomena [ll|, [3, 113, |33| . As 
the turbulence cascades through phase-space, the exci- 
tation of fluctuations at spatial scale £ induces velocity 
structure of scale Sv± in the perpendicular velocity space 
which corresponds to the difference of the Larmor radii 
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FIG. 1. (Color online) Schematic view of the nonlinear 
phase mixing: when the fluctuation scale I is comparable to 
or smaller than the Larmor radius yo, the gyroaverage of the 
electric field induces a decorrelation of the distribution func- 
tion at the velocity-space scale corresponding to the difference 
in Larmor radii — 5v_l/Q. ~ I. 
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FIG 



2. (Color online) Two-dimensional spectral density 
iq\W {k±^p) /Wtot] from one of the forward cascade simu- 



log 

lation reported in 17]. Kinetic turbulence proceeds in the 
position and velocity space simultaneously. 



4 = 5v^/Vt - ^ (see Fig. □ and Refs. [ii-[i3)- In 
other words, when the spatial decorrelation scale is ^, 
two particles with Larmor radii separated by become 
decorrelated, since the gyroaveraged potentials these two 
particles feel are different. This nonlinear phase-space 
mixing_effect was first pointed out by Borland and Ham- 
mett 

In Refs. [l6|, [l^ numerical simulations focused on the 
forward cascade of entropy showed a simultaneous cre- 
ation of structures in position and velocity space in accor- 
dance with the theoretical prediction [ll|, • Figure 
[2] shows a normalized, time-averaged spectral density of 
entropy 
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is characterized by the Hankel transform [l5|, [sj 

g{k,p) = J Jo{pv±)g{k,v±,v\\)dv, (7) 

and "denotes the Fourier-Hankel coefficients. 

In analogy with the Reynolds number in fluid turbu- 
lence, we may introduce an amplitude-dependent dimen- 
sionless number D [16], the ratio of collision time to non- 
linear decorrelation time measured at the thermal Lar- 
mor radius, which characterizes the smallest scales cre- 
ated in both position and velocity space by D~^^^. Thus, 
D quantifies both how "inertial" the turbulence is (in the 
sense of the Reynolds number) and how "kinetic" it is, 
because it measures the nonlinear turnover at the ther- 
mal Larmor radius, which marks the beginning of the 
"nonlinear phase-mixing range." The degree of freedom, 
corresponding to computational problem size, may also 
be characterized by and scales as D^^^ in three phase- 
space dimensions, consisting of two position-space and 
one velocity-space dimensions. 

We note here that the statistical description of our 
phase-space turbulence requires a 2D spectral space 
(/cx,p), in contrast to isotropic fluid turbulence, which 
requires only the scalar wave number. We will refer to 

this 2D spectrum W in the following sections. 



III. THEORY OF FREELY DECAYING 
TURBULENCE 

A. Dual cascade 

With the use of the Hankel transform ([7]) for veloc- 
ity space and the conventional Fourier decomposition for 
position space, we may discuss the evolution of Fourier- 
Hankel modes of g in {k±^p) space. Figure [3] depicts 

the simplest example of the evolution of W{k±^p) in the 
freely decaying turbulence. 

By applying the definition ([7j) to Eqn. ([3|), we find 
that the Fourier-Hankel modes g(k^p) with k^p ^ pvt\^ 
(purple shaded regions marked by "no E"" in Fig. [3]) give 
no contribution to (f){k) due to the orthogonality of the 
Bessel functions. In these regions, therefore, fluctuations 
associated with this part of W can in principle cascade 
to small scales with no effect on the spectral distribution 
oiE. 

On the other hand, the spectrum of E satisfies [HI 
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in a forward cascade simulation, where the v± structure 



where W is given by Eqn. ([6j), i.e., while W is the sum 
over the entire (/c^,p)-plane, the second invariant E is 
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FIG. 3. (Color online) Schematic view of the freely decaying 
turbulence in the {k±,p) space. A diagonal Fourier- Hankel 
mode (indicated by the red circle) cascades towards small 
scales in position and velocity space according to the forward 
cascade of W (blue dotted arrows). Nonlinear interaction 
conserves both W and so the forward cascade must be 
accompanied by the excitation of larger-scale modes on the 
diagonal (green solid arrow) , which corresponds to the inverse 
cascade of E. 



only composed of the "diagonal" k±p = pvn^ compo- 
nents (dashed line in Fig. [3]). In the small-scale limit 
{k±p ^ 1), we may approximate Fq ~ l/{k±p)^ so Eqn. 
dHI) implies 



(9) 



which is indicated for the diagonal components in Fig. [3l 
As we will see, the decaying turbulence evolves to a 
state of local cascade. As a consequence of Eqn. (|9]), 
this state is characterized by an inverse cascade of E and 
forward cascade of W, which can be explained as follows 

The energy of a diagonal Fourier-Hankel mode (indi- 
cated by the red circle) cascades towards small scales in 
position and velocity space according to the forward cas- 
cade of W (blue dotted arrows). Nonlinear interaction 
conserves both W and E, so the forward cascade must 
be accompanied by the excitation of larger-scale modes 
on the diagonal (green solid arrow), which corresponds 
to the inverse cascade of E. Simulation results reported 
in Sec. IIVI indicate that this is the universal property 
of the time-asymptotic state of freely decaying turbu- 
lence. However, there are other interactions possible in 
the transient and/or driven cases (for more details on 
such processes, see Refs. dllsi). 



B. Decay laws 

In this section, we derive scaling laws for the decay of 
collisionless invariants based on some simple phenomeno- 



logical arguments. 

The underlying assumptions are as follows. First, we 
assume that collisionless invariants are dominated by a 
single scale in the same way energy in fluid turbulence 
is dominated by the "energy-containing scale." Then, in 
terms of the amplitude at this scale. 
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where and ip^ are the rms values of the distribution 
function and potential associated with the scale U. Sec- 
ond, we assume that the evolution of is governed by 
the inverse cascade along the diagonal k±p ^ P^th, so 
defining k^ := 1//*, we have 



W - k^pE. 



(11) 



Both assumptions are found to be valid in the numerical 
simulations described in Sec. llVi 

Depending on the strength of collisions compared to 
that of turbulent dynamics, we may derive several differ- 
ent scaling laws. In order to quantify the collisionality, 
we characterize the instantaneous turbulent state via a 
sub-Larmor version of the dimensionless number intro- 
duced in [l^: It is the ratio of the collision time to the 
nonlinear decorrelation time at the scale M, 
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We have taken the collisional decay rate to scale as vklp^ 
because the GK collision operator is second order in ve- 
locity and spatial derivatives [30, 31] and k±p ^ pvth- 
Note that in going from the second to the third expres- 
sion in Eqn. (p!2)) . we used Eqn. (pTj) and 
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valid in the k^p ^ 1 regime from the form of the convec- 
tive derivative {c/ Bq)V {(p) r - V . Note that the factor of 
(/c*p)~^/^ is introduced due to the large argument expan- 
sion of the Bessel function associated with the gyroaver- 
age of (p*. In analogy with the microscopic Reynolds 
number [2^, the initial value and the time evolution of 
provide a natural way to classify various physical 
regimes. 

In the following paragraphs, we describe three differ- 
ent decay laws, classified using Eqn. ([12]). They are the 
weakly collisional (I)* ^ ^o), marginal {D^ ~ Dq) and 
strongly collisional {D^ <C Dq) cases. Here the con- 
stant 1^0 denotes the value of that divide these three 
regimes, which corresponds to the kinetic version of the 
"critical Reynolds number" in Ref. [28]. 

a. Weakly collisional case In the asymptotic limit 
where collision frequency becomes negligible (z^ ^ 0), the 
second invariant E does not decay at all as it is trans- 
ferred to larger scale where dissipation is inactive. On 
the other hand, the first invariant W is transferred to 
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smaller scales, and is dissipated by small but finite col- 
lisions there. The decay rate of W is determined by the 
rate of its transfer to small scales, thus 



dE 
~dt 



0, 



dW 
"dT 



w_ 



(14) 



where we used for the characteristic time of nonlinear 
transfer at scale [see Eqn. (p!3|) ]. Then from ^ t and 

E (X (fl (X , we obtain 



E ~ const.. 



W (X (X t 



-2/3 



(15) 



TABLE I. Index of the runs described in Sec. HVl 



Run z^Tinit kop Nx x Ny x 2Nx init. cond. 



A 


1.7 X 10- 


-3 


15 


64^ 


48^ 


Eqn. 




B 


4.2 X 10" 


-4 


25 


128^ 


96^ 


Eqn. 


([221) 


C 


3.3 X 10" 


-4 


25 


128^ 


96^ 


Eqn. 


(Em 


D 


5.2 X 10" 


-5 


40 


256^ 


192^ 


Eqn. 


(Em 


E 


4.2 X 10" 


-4 


40 


256^ 


192^ 


Eqn. 


(El 


F 


3.8 X 10" 


-4 


25 


256^ 


192^ 


Eqn. 


([241) 



where we used Eqns. (fTTl) and ([T3l) . In this limit Eqn. 
([12]) implies that increases in time as (x t^l'^ . 

b. Marginal case As collision frequency becomes 
large, collisional damping at scale becomes important. 
When it balances with the nonlinear transfer, the two 
terms in 



dW 

~dr 



-—-vklo^W. 



(16) 



become comparable, where we have taken the collisional 
decay rate to scale as fklp^. From Eqns. (fTOj) . ([T3)) and 
^ we obtain the following decay laws of collisionless 
invariants 



E (X k^ (X t 
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W (xt- 
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Substitution of Eqn. ([T7|) into Eqn. ([12]) immediately 
leads to a constant D^. It is noted that Eqn. ([TH]) im- 
plies that the second invariant also decays by collisions 
in a consistent manner: 



dE o o ^ 

dt *^ 



(18) 



c. Strongly collisional case In this case we may re- 
gard the turbulence to be fairly dissipative. We find that 
W and E do not individually satisfy decay laws as powers 
of t. However, Eqn. ([18]) applies, as does the analogous 
equation describing the collisional decay of W ^ 



dW 

~dr 



-vklp^W. 



(19) 



These equations imply [with the help of Eqn. ([TT]) ] decay 
laws for both k^ and the ratio E/W: 



k^ (X t 
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E 
W 



(X 



(20) 



In this case we can deduce that decays in time. It is 
noted that although Eqn. ([20]) holds for both marginal 
and strongly collisional cases, the individual decay laws 
for W and E may vary with [i.e. Eqn. ([TT]) is not 
satisfied here]; however, the decay law for the ratio E/W 
is robust. 



IV. SIMULATION RESULTS 

In this section we show the results of numerical sim- 
ulations performed using the MPI-parallelized nonlinear 
gyrokinetic code AstroGK [s^ • All simulations are made 
in two spatial dimensions {x and y) and two velocity di- 
mensions (energy e = v'^ and pitch angle X = v'^/s) (29| . 



The system size is restricted to = Ly = 27rp, so as to 
focus on the sub-Larmor regime. Time is normalized by 
the initial turnover time 



Tinit 



27r5o 



ckl\\{^init)R\ 



(21) 



where \\{^)n\\ = [{l/uo) JJ \{^)r\^Fo dRdv]'/^ , and 
is the wave number at which initial spectrum is peaked 
[see Eqns. ([22]) and ([24]) below]. Our biggest run (Run 
D in Table H]) used 9,216 processor cores for about 50 
wall-clock hours. 



A. Initial conditions 

In order to investigate the freely decaying turbulence, 
we prepared initial conditions peaked at a high wave 
number ko {kop ^ 1) and made a series of simulations for 
varying initial wave number ko and collision frequency u. 
Six runs were made for two kinds of velocity distribu- 
tion functions (described below) and are listed in Table 
m The different initial velocity distributions allow us to 
vary the initial ratio of invariants. For each initial veloc- 
ity distribution, we made weakly and strongly collisional 
simulations. 

a. Coherent velocity distribution (Runs A-D) The 
first type of the initial velocity distribution is Bessel-like, 
with a Maxwellian envelope Fq, 



g{k,v±,v\\) = ^OT^exp 



k± - ko 



Jo 



k±v± 



Fo. 



(22) 

where the width of the wave- number peak is /cwP = 1- 
From quasi- neutrality [Eqn. ([3])], it can be deduced that 
small-scale velocity oscillation whose period in velocity 
space is comparable to ^'th/(^oP) = ^/^o are needed 
to produce a finite potential (see the discussions in Sec. 
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IIII Ap . Indeed, such oscillatory structure can be found in 
the eigenmodes of the entropy mode [38], which can be 
unstable at quite large k^p. The Fourier-Hankel spec- 
tral density corresponding to Eqn. ([22]) is concentrated 
at k±p = pvth = kop as plotted in the left panel of Fig. 
[4] corresponding to t/rinit = 0, and represented by the 
red spot in Fig. [3l For Eqn. (|22]) the initial ratio of the 
invariants is 



kop, 



(23) 



which we note is the same as what is predicted for the 
time-asymptotic state, given by Eqn. (pTj) . 

b. Random velocity distribution (Runs E and F) 
The second type is a random velocity distribution that 
represents velocity scales 5v/vt\^ ^ l/(A:op), 



^2 



J- ^(2(5, - l)^pJm;MPjV±)Fo, (24) 



TV 

where the width of the wave- number peak is k^p = 1, 
Pj Vth = {k± -\- k^)pr]j^ rjj and 5j are homogeneous ran- 
dom numbers in (0,1), and A/" = 50 is the number of 
random modes for each k^. The factor of ^Jpjn^ is in- 
troduced to cancel the same factor of the Bessel function 
in the asymptotic regime (pj^th ^1). As is shown in 
the left panel (t/rinit = 0) of Fig. El Eqn. corre- 
sponds to a high-density band parallel to the p axis in 
the {k±^p) spectral space. In this case the initial ratio of 
the invariants is 



KqP . 



(25) 



Note that while Eqn. ([22]) is used to represent the veloc- 
ity structure of a coherent mode, the distribution ([24]) 
is designed to mimic the random nature of the velocity 
space that develops from forward cascade (see Ref. jl6i]). 



B. Spectral evolution 

The time evolution of the 2D spectrum W{k±,p) [see 
Eqn. (j6j)] for Runs D and F (Table is shown in Figs.gj 
andO respectively. 

In Fig. m the spectrum is concentrated around k±p = 
pvth = 40 at t = [see also Eqn. (|22|) and Table Hj. The 
spectral density is transferred diagonally to the lower-left 
corner of the {k±^p) space. Since the high-(/c^,p) compo- 
nents suffer strong collisional dissipation, the upper-right 
energy content damps quickly. The remaining lower-left 
transfer dominates after t/rinit ^ 20, and inverse cascade 
follows. 

Nonlinear transfer of the invariants may be directly 
monitored as is done for the forward cascade simulation 



in Refs. [l7|, l39|. Following Ref. [40|, we define a shell 
filtered function by 



keK, 

gK{R) := 9{k)e^'' '', 



(26) 
(27) 



where IC = {k:Kp-l/2< \k\p < Kp + 1/2}. Then the 
evolution of energy in shell K is described as 



d noqf 
dt 2Toi 



■ ^(l-ro)|(^(fe)|2 =^T(^)(if,Q)-collisions, 
fce/c Q 

(28) 

where we introduced an energy transfer function 

T(^HK,Q) := --g^ 1 1 {v>K)R{{^Q)R,9}dRdv, 

(29) 

which measures the rate of energy transferred from shell 
Q to shell K. Here V denotes the spatial volume of the 
domain. The entropy transfer function is defined in a 
similar manner in Ref. jll7|] and is recaptured here in the 
present notation: 



BnV 



gK{{^)R,gQ} 

Fo 



dRdv. (30) 



Note that the shell filtering is performed on (p and g in 
Eqns. dMl) and (|30|), respectively, so that T^^) and T^^) 
both satisfy antisymmetry under exchange of K and Q. 

Snapshots of the normalized energy transfer func- 
tion T^^\K^Q)/E and entropy transfer function 
T^^\K,Q)/W at t/rinit = 77 are shown for the sim- 
ulation of coherent initial condition (Run D) in Fig. [6l 
At this time the peak of the wave-number spectra (not 
shown) is located at /c^p ~ 4. Figure [Ha) shows that the 
energy transfer is (1) well localized along the diagonal 
(meaning local-scale interaction), (2) has strong positive 
values at {Kp.Qp) = (3,4), (2,4), (2,3), (3) has cor- 
responding negative values at {Kp^Qp) = (4,3), (4,2), 
(3,2), and (4) disappears rather quickly at high wave- 
number shells; namely, the energy is transferred from 
large wave-number shells (small scales) to small wave- 
number shells (large scales) around the spectral peak, 
showing clear evidence of the inverse cascade of E. This 
inverse transfer of E creates the peaked, high-density re- 
gion at the diagonal of {k±^p) space, which propagates 
along the diagonal toward the small {k±^p) regime as 
seen at t/rinit = 19 and 96 of Fig. [H At an initial 
transient stage, however, we observe differences including 
nonlocal transfer, which is discussed elsewhere [l8|, [35| . 
Note that as the contribution to E only comes from the 
k±p = pvtii component of the distribution function, the 
transfer in velocity space proceeds in conjunction with 
that of position space, effectively "unwinding" fine struc- 
ture in position and velocity space simultaneously. This 
is a striking feature of the phase-space cascade. 
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FIG. 4. (Color online) Time evolution of the 2D spectra log^Q^ {k_L,p) /W\ for Run D [see Eqn. Q and Table H]. Diagonal 
components (k±p = pvth) are indicated by dotted lines. 
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FIG. 5. (Color online) Time evolution of the 2D spectra \og-^Q[W {k±, p) /W] for Run F [see Eqn. (|6j) and Table [l]. The diagonal 
{k±p = pvth) is indicated by dotted lines. 



On the other hand, from Fig.[6l^b), the entropy transfer 
(1) is well localized along the diagonal (meaning local- 
scale interaction), (2) has a positive peak at {Kp^Qp) = 

(5.4) and corresponding negative peak at {Kp^Qp) = 

(4.5) , (3) extends to larger wave- number shells contrary 
to the transfer of T^^^; namely, the entropy is mostly 
transferred from small wave-number shells (large scales) 
to large wave- number shells (small scales), showing clear 
evidence of the forward cascade of W. This forward 
transfer of W creates the broad off-diagonal spectra seen 
at t/rinit = 19 and 96 of Fig. [H similar to the for- 
ward cascade simulation [17] (see also Fig. [2]). How- 
ever, it is interesting to note the reversed coloring around 
{Kp^Qp) :^ (2,5) of Fig. [6Kb), which is located outside 
of the closest diagonal grids that show forward cascade 
of W. This denotes the inverse transfer of W associated 
with the strong inverse cascade of E; however, its magni- 
tude is less than 1/3 of the peak of the forward transfer 
of W [note also the difference of the scale on the color 



bar in Fig. (6]^ a) and (b)]. 

In the case of random initial condition, the initial ran- 
dom velocity distribution (f24|) is characterized by a verti- 
cal band in {k±,p) space (see Fig. [5]), covering p'^th ^ 25 
(note that kop = 25 as shown in Table [U Run F). Most of 
the energy in the band is transferred to higher wave num- 
ber and then dissipated by strong collisions as it is not 
associated with E (recall the discussions in Sec. IIII A\i . 
For larger scales {k±p < 25), excitation of Fourier-Hankel 
modes is concentrated about the diagonal component 
{k±p = pvth)^ as seen at t/rinit = 13 and 65. As can be 
expected from the time evolution of the 2D spectra (com- 
pare Figs, m and [5]), random-initial-condition cases and 
the coherent-initial-condition cases show similar transfer 
after the transient phase (that is, the transfer functions 
resemble those of Fig. [6|). 

Figure [3 shows several slices of Fig.[5l These slices are 
taken at the peak of the 2D spectra at each time. Ex- 
cept at t/rinit = 65, the diagonal component pvn^ k±p 
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FIG. 6. (Color online) Snapshot of the (a) energy trans- 
fer function T^^\K,Q)/E and (b) entropy transfer function 
T^^\K, Q)/W at t/rinit = 77 for Run D [see Eqns. (|29]), {jSO]) 
and Table II]. 
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FIG. 8. (Color online) Time evolution of for the runs 
indexed in Table U Runs A and E correspond to the strongly 
collisional case, B corresponds to the marginal case and C, D 
and F correspond to the weakly collisional case. 



is excited spontaneously with a nearly Gaussian form, 
which is in contrast to the initial p-spectrum consisting 
of a broad band occupying pvth ^ 25. Thus we con- 
clude that the peaked excitation about the diagonal is 
not merely a reflection of a similarly peaked initial p- 
spectrum. Furthermore, the peaking of the spectrum 
around pvth = k±p^ combined with Eqn. ([9|), justifies the 
approximation given in Eqn. (pT]) . On the other hand, the 
Gaussian peak is surrounded by a fairly broad spectrum 
of an order of magnitude smaller amplitude, which is due 
to the excitation of random velocity fluctuation arising 
from the small-amplitude forward cascade [see Fig.[6l^b)]. 
Note that at t/rinit = 65, the spectrum has a fairly broad 
peak because the scale of the peak is approaching the 
system size (due to a background Maxwellian with the 
thermal velocity '^th)- 

In both cases, the spectra share some qualitative fea- 
tures with the runs reported in Ref. [17] (see also Fig.[2j): 
At large k± and p, the spectral density is broadly dis- 
tributed about the diagonal, which is expected from the 
forward cascade of W However, in each case there 
is a highly peaked diagonal component (whose scale is 
denoted by /* in Sec. IIIIBp that tends to be the domi- 
nant contribution to the total value of the invariants at 
the later stage. This is the component generated by the 
inverse cascade and is discussed in more detail in the 
following sections. 



C. Decay laws 

The time evolution of the dimensionless number 
is shown for each of the runs in Fig. [51 Depending on 
the long-time behavior, we may classify the runs into 
strongly collisional, marginal and weakly collisional cases 
as described in Sec. IIII B[ Run B corresponds to the 
marginal case as approaches a constant (Z^o — 12) 
in t/rinit ^ 30. Runs A and E are strongly colhsional as 
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FIG. 9. (Color online) Time evolution of the ratio of colli- 
sionless invariants for the strongly collisional (Runs A and E) 
and marginal (Run B) cases. The decay law ()20p is drawn for 
comparison. 
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Decay law of each invariant ()17p are drawn for the marginal 
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case. 



decreases in time and Runs C, D and F are weakly 
collisional as increases in time. The evolution of 
differs among weakly collisional runs but approaches the 
theoretical limit oc t^l'^ as increases. 

Figure [9] shows the time evolution of the ratio of two 
collisionless invariants for strongly collisional (Runs A 
and E) and marginal (Run B) cases. Initially, the ratio 
is of the order of l/(A:op) or Xjik^p^) depending on the 
initial velocity distributions [see Eqns. ([23]) and (|25]) ]. 
Coherent cases show an initial phase of constant ratio 
(^/Tinit ^ 6 for Run A and t/rinit ^ 8 for Run B), which 
stems from the fact that the initial condition ([22]) is al- 
most monochromatic with high-(/c^,p), and that both 
collisionless invariants initially decay at the same rate due 
to the collisional damping of the distribution function at 
the scale /c^^. The decay law ([20]) seems consistent with 
both strongly collisional and marginal cases at the later 
stage. The case with the random initial condition tends 
to take a longer time to approach the theoretical line as 
the initial ratio, E jW ~ 1/(/cqP^) [see Eqn. ([25]) ]. is much 
smaller than the asymptotic ratio l/{k^p). A slight de- 
viation from the power law is observed at the last stage 
of the simulation for Runs A and B, which is due to the 
fact that the cascade has reached the largest wave length 
of the system around E /W ^ 0.2, due to the finite size 
of the simulation box. 

The time evolution of individual collisionless invari- 
ants are shown in Fig. [TOjfor the marginal (Run B) and 
weakly collisional (Run D) cases. The marginal run (B) 
shows a reasonable agreement with the theoretical ex- 
pectation ([TT]) . The weakly collisional run shows a sig- 
nificantly slower decay than the marginal case, but still 
faster than the asymptotic {v 0) limit ([T5]) . Compu- 
tational resources beyond those available for this study 
will be needed to unequivocally confirm the asymptotic 
decay laws ([T5]) . 



D. Self-similarity of spectra 

The self-similarity of the spectrum may be investigated 
with the approach of Chasnov for 2D fluid turbulence 
[28j . Chasnov defined an instantaneous length scale of 
the decaying turbulence in terms of the ratio of the two 
invariants, energy (u^) and enstrophy (cj^). This is pos- 
sible because of the relationship that hold between them 
at each scale, namely vorticity is a first-order deriva- 
tive of velocity u. Here we may define an instantaneous 
length scale in terms of the ratio of the two collisionless 
invariants W and E, based on Eqn. (|9|). In general W 
and E can be nearly independent due to the extra degree 
of freedom arising from the velocity space; however, the 
present argument is valid when the spectral density is 
sufficiently concentrated along the diagonal, k±p = pvth- 

We normalize the spectra in velocity space as well as in 
position space. On dimensional grounds, we may define 



Wk{k) = 



W 

K.E{k^,t) 
E 

w ' 



(31) 
(32) 
(33) 



where /c* and are the inverses of the characteristic 
scale length in position and velocity space, respectively, 
determined by [see Eqn. ([TT]) ]. 



w_ 



w 

wth-E" 



(34) 



the wave number k± and velocity wave number p are 
normalized to /c* and p^ , 



k^ 



P 

P = — , 



(35) 
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and the wave-number and velocity-space spectra are de- 
fined by 



y f ToMk )\' 
~f J 2Fo 



\k\=k_ 



(36) 
(37) 



and Eqn. dH). 

When we apply the normalization ([3T]) - ([33|) to simu- 
lation results, one would expect a good coincidence for 
the same value of D^. Namely, if the theory is correct, 
the normalized spectra at different times should collapse 
in the time-asymptotic regime of the marginal case (Run 
B, see also Fig. [8]). 

We first show the normalized spectra defined by Eqns. 
(j3l])-(|33l) in Fig. H] for Run B (see Table IH). As one can 
easily expect from the 2D spectra (Fig. |4]) and the ratio 
(Fig.[9|), the spectral peak moves towards larger scales [to 
smaller {k±^p)] and we may expect self-similar spectra 
in the later stage of the simulation (From Figs. [SHTOl we 
see that the asymptotic regime is attained in the range 
30 ^ t/rinit ^ 150). After the ratio E/W approaches 
the power- law behavior (t/rinit ^ 50, see Fig. [9]), wave- 
number spectra show promising coincidence which indi- 
cates a good self-similarity of the decaying turbulence. 
Notice especially the amazing coincidence of the tail at 
t/rinit = 61 and 120, both in the time- asymptotic regime 
as shown in Figs. [SHTOl We also confirm that the peak 
of the spectra moves towards large scales roughly with 
the scaling oc in this time regime. The last one 

{t/rinit — 240) is offset by a small amount because of 
the fact that the spectral peak has reached the system 
size. The wave- number spectra at the right of the peak 
show slopes somewhat steeper than the theoretical pre- 
diction of the forward cascade [l6|, [l^ • This agrees with 
the expectation from 2D fluid turbulence [28|, as the di- 
mensionless number is not asymptotically large in the 
marginal case. 

As becomes large, the spectral slope becomes shal- 
lower and the forward cascade spectra with the slope 
-4/3 for W and -10/3 for E (see Ref. ^) is expected 
to be recovered. Figure [12] shows the normalized spec- 
tra for Run D. Comparison with Fig. [11] clearly indicates 
that the spectral slopes in Fig. [12] are indeed closer to 
the theoretical expectation of the forward cascade, for 
both spectra of W]^ and E]^. Collapse of the spectra is as 
good as Fig. [TTJ However, as continues to grow for 
V'^init ^ 30 (see Fig. [8]), the spectral slope does change 
slightly, gradually approaching the theoretical prediction 
(W'fc ^ k-^^^ and Ek ^ k'^^^^) as shown in Fig. [12 At 
V^init = 380 the spectral peak has reached the system 
size and the spectra become offset. 
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FIG. 11. (Color online) Normalized spectra for Run B (see 
Table [I|). 



V. SUMMARY 

We presented theoretical and numerical investigations 
of electrostatic, freely decaying turbulence of weakly- 
collisional, magnetized plasmas using the gyrokinetic 
model in 4D phase space (two position-space and two 
velocity-space dimensions). Landau damping was re- 
moved from the system by ignoring variation along the 
background magnetic field. Nonlinear interactions intro- 
duce an amplitude-dependent perpendicular phase mix- 
ing of the gyrophase-independent part of the perturbed 
distribution function, which creates structure in v±_ com- 



11 




0.05 0.1 0.2 0.5 1 2 5 10 20 



k 




0.05 0.1 0.2 0.5 1 2 5 10 20 



k 











—4/3 
































t/ri,it=1 9 

^^init=190 
t/rjnjt=380 











0.05 0.1 0.2 0.5 1 2 5 10 20 

FIG. 12. (Color online) Normalized spectra for Run D (see 
Table m. 

parable in size to spatial structure (see Fig. [1]). 

Since our 2D (in position space) system possesses two 
collisionless invariants [entropy, or free energy, see Eqn. 
(|4]); and energy, see Eqn. (j5j)], a dual cascade (forward 



and inverse cascades) takes place when the initial con- 
dition consisted of small-scale fluctuations in position as 
well as in velocity space such as in Eqns. ([22]) and ([24]) . 
As the dual cascade proceeds, the peak of the spectra 
moves towards large scales in both position and velocity 
space as shown in Figs. IIHSI Nonlinear transfer is diag- 
nosed by the direct numerical simulation, which shows 
a clear evidence of inverse (forward) transfer of energy 
(entropy) (see Fig. [6]). In the inverse cascade, the ve- 
locity space spectrum is highly focused due to the fact 
that energy comes from coherent structure in the velocity 
space [see Eqn. (|9]) and Fig. [7], which is in contrast to 
the broad distribution of velocity scales excited at each 
wave number in the forward cascade [l6|, [l^ • 

Following an example from 2D Navier- Stokes turbu- 
lence a phenomenological theory of decay is pre- 
sented (Sec. nil Bp as well as the numerical simulation 
(Sec. IIV C]) . Several types of asymptotic decay have been 
identified in numerical simulation, which match up well 
with the phenomenological theory using a classification 
based on the kinetic dimensionless number [see Eqn. 
(p!2]) ]. When takes a marginal value, decay laws of 
both invariants are identified [see Eqn. (pT|) . Figs. [8] and 
[To] . In the weakly collisional regime the invariants decay 
more slowly [see Fig. [10]; and in the asymptotic limit 
where the collision frequency becomes negligible (but fi- 
nite), the entropy (or free energy) decays as while 
the energy stays constant [see Eqn. ([T5]) ]. 

In this paper we focused on the time-asymptotic regime 
of the freely decaying turbulence. Although there is a 
range of behavior depending on the strength of dissi- 
pation, the cases are unified by some common features. 
The most prominent of these features is the dual cas- 
cade, whereby the invariant E cascades inversely to large 
scales while the invariant W cascades to small scales. The 
transient and driven cases, hosts to a broader range of 
phenomena, have recently been explored in other works 
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